In [1]:
from __future__ import print_function, division
In [2]:
# This changes the current directory to the base saga directory - make sure to run this first!
# This is necessary to be able to import the py files and use the right directories,
# while keeping all the notebooks in their own directory.
import os
import sys
from time import time
if 'saga_base_dir' not in locals():
saga_base_dir = os.path.abspath('..')
if saga_base_dir not in sys.path:
os.chdir(saga_base_dir)
In [3]:
import hosts
import targeting
import numpy as np
from scipy import interpolate
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy import table
from astropy.table import Table
from astropy.io import fits
from astropy.utils.console import ProgressBar
from collections import Counter
In [4]:
%matplotlib inline
from matplotlib import style, pyplot as plt
plt.style.use('seaborn-deep')
plt.rcParams['image.cmap'] = 'viridis'
plt.rcParams['image.origin'] = 'lower'
plt.rcParams['figure.figsize'] = (14, 8)
plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 16
plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 14
In [5]:
from IPython import display
In [6]:
from decals import make_cutout_comparison_table, fluxivar_to_mag_magerr, compute_sb, DECALS_AP_SIZES, subselect_aperture
In [7]:
hsts = hosts.get_saga_hosts_from_google(clientsecretjsonorfn='client_secrets.json', useobservingsummary=False)
anak = [h for h in hsts if h.name=='AnaK']
assert len(anak)==1
anak = anak[0]
In [8]:
bricknames = []
with open('decals_dr3/anakbricks') as f:
for l in f:
l = l.strip()
if l != '':
bricknames.append(l)
print(bricknames)
In [9]:
base_url = 'http://portal.nersc.gov/project/cosmo/data/legacysurvey/dr3/tractor/{first3}/tractor-{brickname}.fits'
for brickname in ProgressBar(bricknames, ipython_widget=True):
url = base_url.format(brickname=brickname, first3=brickname[:3])
target = os.path.join('decals_dr3/catalogs/', url.split('/')[-1])
if not os.path.isfile(target):
!wget $url -O $target
else:
print(target, 'already exists, not downloading')
In [10]:
bricks = Table.read('decals_dr3/survey-bricks.fits.gz')
bricksdr3 = Table.read('decals_dr3/survey-bricks-dr3.fits.gz')
In [11]:
catalog_fns = ['decals_dr3/catalogs/tractor-{}.fits'.format(bnm) for bnm in bricknames]
decals_catalogs = [Table.read(fn) for fn in catalog_fns]
dcatall = table.vstack(decals_catalogs, metadata_conflicts='silent')
In [12]:
sdss_catalog = Table.read('catalogs/base_sql_nsa{}.fits.gz'.format(anak.nsaid))
In [13]:
plt.scatter(dcatall['ra'], dcatall['dec'], lw=0, c='g', s=2, alpha=.25)
plt.scatter(sdss_catalog['ra'], sdss_catalog['dec'], lw=0, c='r', s=2, alpha=.25)
csize = 1*u.deg
phase = np.arange(360)*u.deg
xc = anak.coords.ra + csize * np.sin(phase)
yc = anak.coords.dec + csize * np.cos(phase)
plt.plot(xc, yc, lw=2, c='k')
Out[13]:
In [14]:
#cut out the non-overlap region
dsc = SkyCoord(dcatall['ra'], dcatall['dec'], unit=u.deg)
dcutall = dcatall[dsc.separation(anak.coords) < 1*u.deg]
plt.scatter(dcutall['ra'], dcutall['dec'], lw=0, c='g', s=2, alpha=.25)
plt.scatter(sdss_catalog['ra'], sdss_catalog['dec'], lw=0, c='r', s=2, alpha=.25)
Out[14]:
Not the over-density in the brick boundaries. I think that's because the bricks have a bit of overlap. Could cut that to get the statistics right if necessary, but for now will leave it in for simplicity.
In [15]:
maggaldepth = -2.5*(np.log10(5*dcutall['decam_galdepth']**-0.5)-9)
maggaldepthr = maggaldepth[:, 2]
In [16]:
plt.figure()
plt.hist(maggaldepthr[np.isfinite(maggaldepthr)], bins=100, histtype='step', log=False)
plt.xlabel('galdepth in $r$')
plt.ylabel('$n$')
plt.tight_layout()
plt.xlim(22, 25.5)
plt.figure()
plt.hist(maggaldepthr[np.isfinite(maggaldepthr)], bins=100, histtype='step', log=True)
plt.xlabel('galdepth in $r$')
plt.ylabel('$n$')
plt.tight_layout()
plt.xlim(22, 25.5)
Out[16]:
In [17]:
r, r_err = fluxivar_to_mag_magerr(dcutall['decam_flux'][:, 2], dcutall['decam_flux_ivar'][:, 2])
In [18]:
plt.hexbin(r, r_err, extent=(15,26,0,1), bins='log',gridsize=200)
Out[18]:
In [19]:
H, xe, ye = np.histogram2d(r, r_err, bins=200, range=((15,26),(0,1)))
lH = np.log(H)
lH[~np.isfinite(lH)] = np.min(lH[np.isfinite(lH)])-1
plt.pcolor(xe, ye, lH.T)
plt.xlim(15, 26)
Out[19]:
In [20]:
for dcat in [dcutall]:
for magnm, idx in zip('grz', [1, 2, 4]):
mag, mag_err = fluxivar_to_mag_magerr(dcat['decam_flux'][:, idx], dcat['decam_flux_ivar'][:, idx])
dcat[magnm] = mag
dcat[magnm + '_err'] = mag_err
dcat['sb_r_0.5'] = compute_sb(0.5*u.arcsec, dcat['decam_apflux'][:, 2, :])
dcat['sb_r_0.75'] = compute_sb(0.75*u.arcsec, dcat['decam_apflux'][:, 2, :])
dcat['sb_r_1'] = compute_sb(1.0*u.arcsec, dcat['decam_apflux'][:, 2, :])
dcat['sb_r_2'] = compute_sb(2.0*u.arcsec, dcat['decam_apflux'][:, 2, :])
Details outlined in "DECALS low-SB_completeness figures" notebook
In [21]:
dcat = dcutall
rap_fluxes = dcat['decam_apflux'][:, 2, :]
expflux_r = np.empty_like(rap_fluxes[:, 0])
rad = np.empty(len(rap_fluxes[:, 0]))
ap_sizesv = DECALS_AP_SIZES.to(u.arcsec).value
intr = interpolate.BarycentricInterpolator(ap_sizesv, [0]*len(ap_sizesv))
for i in ProgressBar(range(len(rap_fluxes)), ipython_widget=True):
f = rap_fluxes[i]
if dcat['type'][i] == 'PSF ':
r = dcat['decam_psfsize'][i, 2]
elif dcat['type'][i] == 'DEV ':
r = dcat['shapeDev_r'][i]
else:
r = dcat['shapeExp_r'][i]
intr.set_yi(f)
expflux_r[i] = intr(r)
rad[i] = r
dcat['sdss_like_sb_r'] = compute_sb(rad*u.arcsec, np.array(expflux_r))
In [22]:
dsc = SkyCoord(dcutall['ra'], dcutall['dec'], unit=u.deg)
ssc = SkyCoord(sdss_catalog['ra'], sdss_catalog['dec'], unit=u.deg)
threshold = 1*u.arcsec
In [23]:
idx, d2d, _ = ssc.match_to_catalog_sky(dsc)
plt.hist(d2d.arcsec, bins=100, range=(0, 3),histtype='step', log=True)
plt.axvline(threshold.to(u.arcsec).value, c='k')
None
In [24]:
dmatchmsk = idx[d2d<threshold]
dmatch = dcutall[dmatchmsk]
smatch = sdss_catalog[d2d<threshold]
In [25]:
idx, d2d, _ = dsc.match_to_catalog_sky(ssc)
dnomatchmsk = d2d>threshold
dnomatch = dcutall[dnomatchmsk]
In [26]:
idx, d2d, _ = ssc.match_to_catalog_sky(dsc[dnomatchmsk])
plt.hist(d2d.arcsec, bins=100, range=(0, 3),histtype='step', log=True)
plt.axvline(1.0, c='k')
None
In [27]:
idx, d2d, _ = dsc[dmatchmsk].match_to_catalog_sky(ssc)
plt.hist(d2d.arcsec, bins=100, range=(0, 3),histtype='step', log=True)
plt.axvline(1.0, c='k')
None
In [28]:
plt.figure(figsize=(12, 10))
xnm = 'r'
ynm = 'sb_r_0.5'
dmextinction = -2.5*np.log10(dmatch['decam_mw_transmission'][:, 2])
x = dmatch[xnm] - dmextinction
y = dmatch[ynm] - dmextinction
dstar = dmatch['type']=='PSF '
sstar = smatch['phot_sg']==6
plt.scatter(x[~dstar&~sstar], y[~dstar&~sstar], c='b', lw=0, alpha=.4, s=2, label='Glx in Both')
plt.scatter(x[dstar&~sstar], y[dstar&~sstar], c='r', lw=0, alpha=.5, s=3, label='Star in DECaLS, Glx in SDSS')
plt.scatter(x[~dstar&sstar], y[~dstar&sstar], c='g', lw=0, alpha=.75, s=5, label='Glx in DECaLS, Star in SDSS')
dnstar = dnomatch['type']=='PSF '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
plt.scatter((dnomatch[xnm] - dnoext)[~dnstar], (dnomatch[ynm] - dnoext)[~dnstar],
c='k', lw=0, alpha=.4, s=2, label='Glx in DECALS, not in SDSS')
#plt.hlines(24.5, 17, 21, colors='k', linestyles='--')
plt.axvline(20.75, color='k', ls=':')
plt.xlim(17, 24.5)
plt.ylim(18, 28)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$SB_{0.5^{\prime \prime}, {\rm DECaLS}}$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc='lower right', fontsize=20)
#plt.savefig(os.environ['HOME'] + '/Desktop/fig_decals.pdf')
Out[28]:
In [29]:
to_check = dcutall[(dcutall['r']<21)&(dcutall['sb_r_0.5']>24.5) & (dcutall['type']!='PSF ')]
np.random.shuffle(to_check)
print('n=', len(to_check))
make_cutout_comparison_table(to_check[:25])
Out[29]:
In [30]:
to_check = dnomatch[((18<dnomatch['r']) & (dnomatch['r']<21) &
(18<dnomatch['sb_r_0.5']) & (dnomatch['sb_r_0.5']<24.5) &
(dnomatch['type']!='PSF '))]
#np.random.shuffle(to_check)
print('n=', len(to_check))
make_cutout_comparison_table(to_check[:25])
Out[30]:
In [31]:
msk = np.zeros(len(dnomatch), dtype=bool)
for nm in ['327464_6040', '327465_5904']:
bid, oid = nm.split('_')
msk |= ((dnomatch['objid']==int(oid))&(dnomatch['brickid']==int(bid)))
make_cutout_comparison_table(dnomatch[msk])
Out[31]:
In [32]:
dnomatch[msk]
Out[32]:
Both have anomalously large sizes... clearly they are not actually r<21
In [33]:
to_check = dmatch[~dstar&sstar&(dmatch['r']<21)]
np.random.shuffle(to_check)
print(len(to_check))
make_cutout_comparison_table(to_check[:25])
Out[33]:
Important thing to understand here: this is not the same as the SDSS SB's for a subtle reason. The SDSS SB's are all from assuming an exponential profile. Here, instead (because DECaLS doesn't have the data for exp profiles if the object is not found to be exp type), we use an SB within either the exponential half-light radius, the de Vaucouleurs effective radius, or the seeing FWHM for exp galaxies, DeV galaxies, or PSF's, respectively.
The flux is also not the modeled flux, but rather an interpolated flux at that radius from the apertures.
In [34]:
plt.figure(figsize=(12, 10))
xnm = 'r'
ynm = 'sdss_like_sb_r'
dmextinction = -2.5*np.log10(dmatch['decam_mw_transmission'][:, 2])
x = dmatch[xnm] - dmextinction
y = dmatch[ynm] - dmextinction
dstar = dmatch['type']=='PSF '
sstar = smatch['phot_sg']==6
dnstar = dnomatch['type']=='PSF '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
plt.scatter((dnomatch[xnm] - dnoext)[~dnstar], (dnomatch[ynm] - dnoext)[~dnstar],
c='k', lw=0, alpha=.25, label='Glx in DECALS, not in SDSS', marker='s')
plt.scatter((dnomatch[xnm] - dnoext)[dnstar], (dnomatch[ynm] - dnoext)[dnstar],
c='k', lw=0, alpha=.15, s=60,
label='Star in DECALS, not in SDSS', marker='*')
plt.scatter(x[dstar&sstar], y[dstar&sstar], c='m', lw=0, alpha=.25, s=60,
label='Star in both', marker='*')
plt.scatter(x[~dstar&~sstar], y[~dstar&~sstar], c='b', lw=0, alpha=.25, label='Glx in Both')
plt.scatter(x[dstar&~sstar], y[dstar&~sstar], c='r', lw=0, alpha=.25,
label='Star in DECaLS, Glx in SDSS', marker='^', s=50)
plt.scatter(x[~dstar&sstar], y[~dstar&sstar], c='g', lw=0, alpha=.25,
label='Glx in DECaLS, Star in SDSS', marker='^', s=50)
plt.axvline(20.75, color='k', ls=':')
plt.xlim(17, 24.5)
plt.ylim(18, 30)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$SB_{sdssish}, {\rm DECaLS}}$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc='lower right', fontsize=20)
#change to "Should have only the Glx both, Glx D/not in S, and Gx D/star in S points"
#OUTPUT THESE DATA FOR MARLA,
Out[34]:
Need to calculate the "dmag"'s - see "DECALS low-SB_completeness Residuals" notebook for more info on this
In [56]:
def extract_dmags(dcat):
apmag, apmagerr = fluxivar_to_mag_magerr(dcat['decam_apflux'], dcat['decam_apflux_ivar'])
apmagres, _ = fluxivar_to_mag_magerr(dcat['decam_apflux_resid'], dcat['decam_apflux_ivar'])
dmags = {}
for band, ap in [('r', 1.0*u.arcsec), ('g', 1.0*u.arcsec), ('r', 5.0*u.arcsec), ('g', 5.0*u.arcsec)]:
nm = 'dmagap_{}_{}asec'.format(band, ap.value)
dmags[nm] = subselect_aperture(apmag, band, ap) - subselect_aperture(apmagres, band, ap)
return dmags
In [59]:
to_send_to_marla = dmatch['ra', 'dec']
to_send_to_marla['r0'] = x
to_send_to_marla['r_err'] = dmatch['r_err']
to_send_to_marla['Rh_exp'] = dmatch['shapeExp_r']
to_send_to_marla['Rh_dev'] = dmatch['shapeDev_r']
to_send_to_marla['Rh_psf'] = dmatch['decam_psfsize']/2.
to_send_to_marla['sb_sdsslike'] = y
to_send_to_marla['id'] = ['{}_{}'.format(row['brickid'], row['objid']) for row in dmatch]
col = to_send_to_marla['id']
del to_send_to_marla['id']
to_send_to_marla.add_column(col, 0)
to_send_to_marla
to_send_to_marla['insdss'] = True
to_send_to_marla['star_decals'] = dstar
to_send_to_marla['star_sdss'] = sstar
for nm, arr in extract_dmags(dmatch).items():
to_send_to_marla[nm] = arr
In [60]:
toadd = dnomatch['ra', 'dec']
toadd['r0'] = dnomatch[xnm] - dnoext
toadd['r_err'] = dnomatch['r_err']
toadd['Rh_exp'] = dnomatch['shapeExp_r']
toadd['Rh_dev'] = dnomatch['shapeDev_r']
toadd['Rh_psf'] = dnomatch['decam_psfsize']/2.
toadd['sb_sdsslike'] = dnomatch[ynm] - dnoext
toadd.add_column(table.Column(name='id', data=['{}_{}'.format(row['brickid'], row['objid']) for row in dnomatch]),
index=0)
toadd['insdss'] = False
toadd['star_decals'] = dnstar
toadd['star_sdss'] = False
for nm, arr in extract_dmags(dnomatch).items():
toadd[nm] = arr
to_send_to_marla = table.vstack([to_send_to_marla, toadd])
In [61]:
to_send_to_marla
Out[61]:
In [63]:
to_send_to_marla.write('/Users/erik/Dropbox/SAGA/paper1_rvssb_data.fits.gz')
The below plot is to see if there's significant differences between the two galaxy types, since their SB calculations are quite different
In [38]:
plt.figure(figsize=(12, 10))
xnm = 'r'
ynm = 'sdss_like_sb_r'
dmextinction = -2.5*np.log10(dmatch['decam_mw_transmission'][:, 2])
x = dmatch[xnm] - dmextinction
y = dmatch[ynm] - dmextinction
dstar = dmatch['type']=='PSF '
ddev = dmatch['type']=='DEV '
sstar = smatch['phot_sg']==6
dnstar = dnomatch['type']=='PSF '
dndev = dnomatch['type']=='DEV '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
plt.scatter((dnomatch[xnm] - dnoext)[~dnstar&~dndev], (dnomatch[ynm] - dnoext)[~dnstar&~dndev],
c='g', lw=0, alpha=.15, label='Exp Glx in DECALS, not in SDSS', marker='*', s=40)
plt.scatter((dnomatch[xnm] - dnoext)[~dnstar&dndev], (dnomatch[ynm] - dnoext)[~dnstar&dndev],
c='k', lw=0, alpha=.15, label='DeV Glx in DECALS, not in SDSS', marker='^', s=30)
plt.scatter(x[~dstar&ddev], y[~dstar&ddev], c='b', lw=0, alpha=.15,
label='Exp Glx in DECaLS, in SDSS', marker='o', s=15)
plt.scatter(x[~dstar&~ddev], y[~dstar&~ddev], c='r', lw=0, alpha=.15,
label='DeV Glx in DECaLS, in SDSS', marker='s', s=20)
#plt.hlines(24.5, 17, 21, colors='k', linestyles='--')
plt.axvline(20.75, color='k', ls=':')
plt.xlim(17, 24.5)
plt.ylim(18, 30)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$SB_{sdssish}, {\rm DECaLS}}$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc='lower right', fontsize=20)
Out[38]:
In [39]:
plt.figure(figsize=(12, 10))
xnm = 'r'
ynm = 'sdss_like_sb_r'
dmextinction = -2.5*np.log10(dmatch['decam_mw_transmission'][:, 2])
x = dmatch[xnm] - dmextinction
y = dmatch[ynm] - dmextinction
dstar = dmatch['type']=='PSF '
sstar = smatch['phot_sg']==6
dnstar = dnomatch['type']=='PSF '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
plt.scatter((dnomatch[xnm] - dnoext)[~dnstar], (dnomatch[ynm] - dnoext)[~dnstar],
c='k', lw=0, alpha=.25, label='Glx in DECALS, not in SDSS', marker='s', s=10)
# plt.scatter((dnomatch[xnm] - dnoext)[dnstar], (dnomatch[ynm] - dnoext)[dnstar],
# c='k', lw=0, alpha=.15, s=60,
# label='Star in DECALS, not in SDSS', marker='*')
# plt.scatter(x[dstar&sstar], y[dstar&sstar], c='m', lw=0, alpha=.25, s=60,
# label='Star in both', marker='*')
plt.scatter(x[~dstar&~sstar], y[~dstar&~sstar], c='b', lw=0, alpha=.25, label='Glx in Both', s=5)
# plt.scatter(x[dstar&~sstar], y[dstar&~sstar], c='r', lw=0, alpha=.25,
# label='Star in DECaLS, Glx in SDSS', marker='^', s=50)
plt.scatter(x[~dstar&sstar], y[~dstar&sstar], c='r', lw=0, alpha=.25,
label='Glx in DECaLS, Star in SDSS', marker='*', s=60)
plt.axvline(20.75, color='k', ls=':')
plt.xlim(17, 24.5)
plt.ylim(18, 30)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$SB_{sdssish}, {\rm DECaLS}}$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc='lower right', fontsize=20)
#plt.savefig(os.environ['HOME'] + '/Desktop/fig_decals.pdf')
Out[39]:
In [40]:
plt.figure(figsize=(12, 10))
xnm = 'r'
dmextinction = -2.5*np.log10(dmatch['decam_mw_transmission'][:, 2])
dstar = dmatch['type']=='PSF '
sstar = smatch['phot_sg']==6
dnstar = dnomatch['type']=='PSF '
dnoext = -2.5*np.log10(dnomatch['decam_mw_transmission'][:, 2])
r_match = dmatch[xnm] - dmextinction
r_nomatch = dnomatch[xnm] - dnoext
r_both = np.concatenate((r_match, r_nomatch))
bothstar = np.concatenate((dstar, dnstar))
plt.hist(r_match[~dstar], histtype='step', color='b', range=(17, 25), bins=100, log=True)
plt.hist(r_both[~bothstar], histtype='step', color='k', range=(17, 25), bins=100, log=True)
plt.xlabel(r'$r_{0, {\rm DECaLS}}$', fontsize=28)
plt.ylabel(r'$n$', fontsize=28)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.tight_layout()
In [96]:
sr = smatch['r']
dr = dmatch['r']
plt.scatter(sr, dr-sr)
plt.ylim(-1,2)
Out[96]: